rm(list = ls())

##Load Packages

library(Amelia)
library(stargazer)
library(VGAM)
library(haven)
library(data.table)
library(estimatr)
library(betareg)
library(mice)
library(texreg)
library(plm)
library(multiwayvcov)
library(lmtest)
library(AER)
library(DataCombine)

#Set Working Directory

#setwd("SET YOUR WORKING DIRECTORY TO DATA")

##Load in Data

fgmc <- read_dta("fgmc.dta")

fgmc$year <- as.factor(fgmc$BirthYearAge)
fgmc$country <- as.factor(fgmc$CountryName)
fgmc$democracy <- ifelse(fgmc$polity2 > 5, 1, 0)
fgmc$autocracy <- ifelse(fgmc$polity2 < -5, 1, 0)

#Convert proportion to percentage
#fgmc$BirthYearAgeFGMCwtd <- (fgmc$BirthYearAgeFGMCwtd * 100)
fgmc$BirthYearAgeFGMCwtdse <- (fgmc$BirthYearAgeFGMCwtdse * 100)

#De-mean and normalize non-binary variables
fgmc$GDPpc <- (fgmc$GDPpc - mean(fgmc$GDPpc, na.rm = T)) / sd(fgmc$GDPpc, na.rm = T)
fgmc$popdensity <- (fgmc$popdensity - mean(fgmc$popdensity, na.rm = T)) / sd(fgmc$popdensity, na.rm = T)
fgmc$durable <- (fgmc$durable - mean(fgmc$durable, na.rm = T)) / sd(fgmc$durable, na.rm = T)
fgmc$enrollmentgenderparityindex <- (fgmc$enrollmentgenderparityindex - mean(fgmc$enrollmentgenderparityindex, na.rm = T)) / sd(fgmc$enrollmentgenderparityindex, na.rm = T)
fgmc$netodagdp <- (fgmc$netodagdp - mean(fgmc$netodagdp, na.rm = T)) / sd(fgmc$netodagdp, na.rm = T)
fgmc$Trade <- (fgmc$Trade - mean(fgmc$Trade, na.rm = T)) / sd(fgmc$Trade, na.rm = T)

#Create Lag/3 Year Average
fgmc <- fgmc[order(fgmc$CountryName, fgmc$BirthYearAge),]

fgmc <- slide(fgmc, Var = "BirthYearAge", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -1, NewVar = "byagelag")
fgmc <- slide(fgmc, Var = "BirthYearAge", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -2, NewVar = "byagelag2")
fgmc <- slide(fgmc, Var = "BirthYearAge", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -3, NewVar = "byagelag3")

fgmc <- slide(fgmc, Var = "BirthYearAgeFGMCwtd", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -1, NewVar = "fgmclag")
fgmc <- slide(fgmc, Var = "BirthYearAgeFGMCwtd", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -2, NewVar = "lag2")
fgmc <- slide(fgmc, Var = "BirthYearAgeFGMCwtd", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -3, NewVar = "lag3")

fgmc <- slide(fgmc, Var = "BirthYearAgeFGMCobs", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -1, NewVar = "obslag")
fgmc <- slide(fgmc, Var = "BirthYearAgeFGMCobs", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -2, NewVar = "obslag2")
fgmc <- slide(fgmc, Var = "BirthYearAgeFGMCobs", GroupVar = "CountryName", TimeVar = "BirthYearAge",
              slideBy = -3, NewVar = "obslag3")

#unweighted
#fgmc$threeyravg <- (fgmc$fgmclag + fgmc$lag2 + fgmc$lag3) / 3

#weighted
fgmc$threeyravg <- (fgmc$fgmclag*fgmc$obslag + fgmc$lag2*fgmc$obslag2 + fgmc$lag3*fgmc$obslag3) / (fgmc$obslag + fgmc$obslag2 + fgmc$obslag3)

#fix non-subsequent years
fgmc$threeyravg <- ifelse((fgmc$byagelag == fgmc$BirthYearAge - 1) & (fgmc$byagelag2 == fgmc$BirthYearAge - 2) & (fgmc$byagelag3 == fgmc$BirthYearAge - 3), fgmc$threeyravg, NA)

#Create Panel Datasets
fgmcp <- pdata.frame(fgmc, index=c("CountryName","BirthYearAge"), drop.index=TRUE, row.names=TRUE)
fgmc85 <- subset(fgmc, BirthYearAge <= 1985)
fgmc85p <- pdata.frame(fgmc85, index=c("CountryName","BirthYearAge"), drop.index=TRUE, row.names=TRUE)

####-------- Cohort Level Unit

## ------ OLS - All data

mod.fe1 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + GDPpc, data = fgmcp, model = "within")
mod.fe2 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + popdensity, data = fgmcp, model = "within")
mod.fe3 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + netodagdp, data = fgmcp, model = "within")
mod.fe4 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + Trade, data = fgmcp, model = "within")
mod.fe5 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + autocracy, data = fgmcp, model = "within")
mod.fe6 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + durable, data = fgmcp, model = "within")
mod.fe7 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + instab, data = fgmcp, model = "within")
mod.fe8 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + war, data = fgmcp, model = "within")
mod.fe9 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + outlaw, data = fgmcp, model = "within")
mod.fe10 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + enrollmentgenderparityindex, data = fgmcp, model = "within")

coeftest(mod.fe1, vcov=vcovHC(mod.fe1,type="HC1",cluster="group"))
coeftest(mod.fe2, vcov=vcovHC(mod.fe2,type="HC1",cluster="group"))
coeftest(mod.fe3, vcov=vcovHC(mod.fe3,type="HC1",cluster="group"))
coeftest(mod.fe4, vcov=vcovHC(mod.fe4,type="HC1",cluster="group"))
coeftest(mod.fe5, vcov=vcovHC(mod.fe5,type="HC1",cluster="group"))
coeftest(mod.fe6, vcov=vcovHC(mod.fe6,type="HC1",cluster="group"))
coeftest(mod.fe7, vcov=vcovHC(mod.fe7,type="HC1",cluster="group"))
coeftest(mod.fe8, vcov=vcovHC(mod.fe8,type="HC1",cluster="group"))
coeftest(mod.fe9, vcov=vcovHC(mod.fe9,type="HC1",cluster="group"))
coeftest(mod.fe10, vcov=vcovHC(mod.fe10,type="HC1",cluster="group"))

htmlreg(list(mod.fe1, mod.fe2, mod.fe3, mod.fe4, mod.fe10), file = "table1",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmclag",
        custom.coef.names = c("GDP Per Capita", "Population Density", "Net ODA", "Trade",
                              "Female Education"),
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

htmlreg(list(mod.fe5, mod.fe6, mod.fe7, mod.fe8, mod.fe9), file = "table2",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmclag",
        custom.coef.names = c("Autocracy", "Regime Durability", "Political Instability",
                              "Civil Conflict", "FGM/C Outlawed"),
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

####### ---------- Repeat with only complete cases

#remove unnecessary variables
keep <- c("BirthYearAgeFGMCwtd", "fgmclag", "GDPpc", "popdensity", "autocracy", "durable",
          "outlaw", "enrollmentgenderparityindex", "netodagdp", "Trade", "instab", "CountryName", "BirthYearAge",
          "threeyravg", "war")
fgmc2 <- fgmc[keep]

fgmc2 <- fgmc2[which(
  is.na(fgmc2$GDPpc) == FALSE &
    is.na(fgmc2$popdensity) == FALSE &
    is.na(fgmc2$autocracy) == FALSE &
    is.na(fgmc2$durable) == FALSE &
    is.na(fgmc2$outlaw) == FALSE &
    is.na(fgmc2$enrollmentgenderparityindex) == FALSE &
    is.na(fgmc2$netodagdp) == FALSE &
    is.na(fgmc2$Trade) == FALSE &
    is.na(fgmc2$instab) == FALSE
),]

fgmc2p <- pdata.frame(fgmc2, index=c("CountryName","BirthYearAge"), drop.index=TRUE, row.names=TRUE)
fgmc285 <- subset(fgmc2, BirthYearAge <= 1985)
fgmc285p <- pdata.frame(fgmc285, index=c("CountryName","BirthYearAge"), drop.index=TRUE, row.names=TRUE)

mod.fe1 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + GDPpc, data = fgmc2p, model = "within")
mod.fe2 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + popdensity, data = fgmc2p, model = "within")
mod.fe3 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + netodagdp, data = fgmc2p, model = "within")
mod.fe4 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + Trade, data = fgmc2p, model = "within")
mod.fe5 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + autocracy, data = fgmc2p, model = "within")
mod.fe6 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + durable, data = fgmc2p, model = "within")
mod.fe7 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + instab, data = fgmc2p, model = "within")
mod.fe8 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + war, data = fgmc2p, model = "within")
mod.fe9 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + outlaw, data = fgmc2p, model = "within")
mod.fe10 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + enrollmentgenderparityindex, data = fgmc2p, model = "within")

htmlreg(list(mod.fe1, mod.fe2, mod.fe3, mod.fe4, mod.fe10), file = "table11.doc",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmclag",
        custom.coef.names = c("GDP Per Capita", "Population Density", "Net ODA", "Trade",
                              "Female Education"),
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

htmlreg(list(mod.fe5, mod.fe6, mod.fe7, mod.fe8, mod.fe9), file = "table12.doc",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmclag",
        custom.coef.names = c("Autocracy", "Regime Durability", "Political Instability",
                              "Civil Conflict", "FGM/C Outlawed"),
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)


coeftest(mod.fe1, vcov=vcovHC(mod.fe1,type="HC1",cluster="group"))
coeftest(mod.fe2, vcov=vcovHC(mod.fe2,type="HC1",cluster="group"))
coeftest(mod.fe3, vcov=vcovHC(mod.fe3,type="HC1",cluster="group"))
coeftest(mod.fe4, vcov=vcovHC(mod.fe4,type="HC1",cluster="group"))
coeftest(mod.fe5, vcov=vcovHC(mod.fe5,type="HC1",cluster="group"))
coeftest(mod.fe6, vcov=vcovHC(mod.fe6,type="HC1",cluster="group"))
coeftest(mod.fe7, vcov=vcovHC(mod.fe7,type="HC1",cluster="group"))
coeftest(mod.fe8, vcov=vcovHC(mod.fe8,type="HC1",cluster="group"))
coeftest(mod.fe9, vcov=vcovHC(mod.fe9,type="HC1",cluster="group"))
coeftest(mod.fe10, vcov=vcovHC(mod.fe10,type="HC1",cluster="group"))

#3 yr lag instead

mod.avg1 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + GDPpc, data = fgmc2p, model = "within")
mod.avg2 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + popdensity, data = fgmc2p, model = "within")
mod.avg3 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + netodagdp, data = fgmc2p, model = "within")
mod.avg4 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + Trade, data = fgmc2p, model = "within")
mod.avg5 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + autocracy, data = fgmc2p, model = "within")
mod.avg6 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + durable, data = fgmc2p, model = "within")
mod.avg7 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + instab, data = fgmc2p, model = "within")
mod.avg8 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + war, data = fgmc2p, model = "within")
mod.avg9 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + outlaw, data = fgmc2p, model = "within")
mod.avg10 <- plm(BirthYearAgeFGMCwtd ~ threeyravg + enrollmentgenderparityindex, data = fgmc2p, model = "within")

htmlreg(list(mod.avg1, mod.avg2, mod.avg3, mod.avg4, mod.avg10), file = "table13.doc",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmc3yravg",
        custom.coef.names = c("GDP Per Capita", "Population Density", "Net ODA", "Trade",
                              "Female Education"),
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

htmlreg(list(mod.avg5, mod.avg6, mod.avg7, mod.avg8, mod.avg9), file = "table14.doc",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmc3yravg",
        custom.coef.names = c("Autocracy", "Regime Durability", "Political Instability",
                              "Civil Conflict", "FGM/C Outlawed"),
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

coeftest(mod.avg1, vcov=vcovHC(mod.avg1,type="HC1",cluster="group"))
coeftest(mod.avg2, vcov=vcovHC(mod.avg2,type="HC1",cluster="group"))
coeftest(mod.avg3, vcov=vcovHC(mod.avg3,type="HC1",cluster="group"))
coeftest(mod.avg4, vcov=vcovHC(mod.avg4,type="HC1",cluster="group"))
coeftest(mod.avg5, vcov=vcovHC(mod.avg5,type="HC1",cluster="group"))
coeftest(mod.avg6, vcov=vcovHC(mod.avg6,type="HC1",cluster="group"))
coeftest(mod.avg7, vcov=vcovHC(mod.avg7,type="HC1",cluster="group"))
coeftest(mod.avg8, vcov=vcovHC(mod.avg8,type="HC1",cluster="group"))
coeftest(mod.avg9, vcov=vcovHC(mod.avg9,type="HC1",cluster="group"))
coeftest(mod.avg10, vcov=vcovHC(mod.avg10,type="HC1",cluster="group"))

#Change Years to eliminate potential programmatic effects

mod.pre1 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + GDPpc, data = fgmc285p, model = "within")
mod.pre2 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + popdensity, data = fgmc285p, model = "within")
mod.pre3 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + netodagdp, data = fgmc285p, model = "within")
mod.pre4 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + Trade, data = fgmc285p, model = "within")
mod.pre5 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + autocracy, data = fgmc285p, model = "within")
mod.pre6 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + durable, data = fgmc285p, model = "within")
mod.pre7 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + instab, data = fgmc285p, model = "within")
mod.pre8 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + war, data = fgmc285p, model = "within")
#mod.pre9 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + outlaw, data = fgmc285p, model = "within")
mod.pre10 <- plm(BirthYearAgeFGMCwtd ~ fgmclag + enrollmentgenderparityindex, data = fgmc285p, model = "within")

htmlreg(list(mod.pre1, mod.pre2, mod.pre3, mod.pre4, mod.pre10), file = "table15.doc",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmclag",
        custom.coef.names = c("GDP Per Capita", "Population Density", "Net ODA", "Trade",
                              "Female Education"),
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

htmlreg(list(mod.pre5, mod.pre6, mod.pre7, mod.pre8), file = "table16.doc",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmclag",
        custom.coef.names = c("Autocracy", "Regime Durability", "Political Instability",
                              "Civil Conflict"),
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

coeftest(mod.pre1, vcov=vcovHC(mod.pre1,type="HC1",cluster="group"))
coeftest(mod.pre2, vcov=vcovHC(mod.pre2,type="HC1",cluster="group"))
coeftest(mod.pre3, vcov=vcovHC(mod.pre3,type="HC1",cluster="group"))
coeftest(mod.pre4, vcov=vcovHC(mod.pre4,type="HC1",cluster="group"))
coeftest(mod.pre5, vcov=vcovHC(mod.pre5,type="HC1",cluster="group"))
coeftest(mod.pre6, vcov=vcovHC(mod.pre6,type="HC1",cluster="group"))
coeftest(mod.pre7, vcov=vcovHC(mod.pre7,type="HC1",cluster="group"))
coeftest(mod.pre8, vcov=vcovHC(mod.pre8,type="HC1",cluster="group"))
coeftest(mod.pre10, vcov=vcovHC(mod.pre10,type="HC1",cluster="group"))

## ------ Beta

fgmc2$BirthYearAgeFGMCwtd2 <- (fgmc2$BirthYearAgeFGMCwtd/100)
bet1 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + GDPpc + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet2 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + popdensity + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet3 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + netodagdp + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet4 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + Trade + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet5 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + enrollmentgenderparityindex  + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet6 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + autocracy + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet7 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + durable + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet8 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + instab + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet9 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + war + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))
bet10 <- betareg(BirthYearAgeFGMCwtd2 ~ fgmclag + outlaw + as.factor(CountryName), data=subset(fgmc2, BirthYearAgeFGMCwtd2 > 0))

htmlreg(list(bet1, bet2, bet3, bet4, bet5), file = "table17.doc",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmclag",
        omit = "as.factor",
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

htmlreg(list(bet6, bet7, bet8, bet9, bet10), file = "table18.doc",
        stars = c(0.01, 0.05, 0.1), omit.coef = "fgmclag",
        custom.note = "Standard Errors Clustered by Country in Parentheses", include.ci = F)

### ----- Create MI Datasets --- WARNING, THIS WILL TAKE A WHILE TO RUN

FGMImpute <- amelia(x=fgmc, m=20, p2s=1,  ts = "BirthYearAge", cs = "CountryName",
                    idvars = c("region_n", "paris", "initprev", "initprevsquared", "fgmclag", "threeyravg",
                               "yearsind", "ccode", "outlaw", "year", "country"),
                    polytime = 1, intercs = F, empri = .1*nrow(fgmc))

#Run models on imputed data

miresults <- function(formula) {
  b.out<-NULL
  se.out<-NULL
  cse.out<-NULL
  lm.out <- NULL
  for(i in 1:FGMImpute$m) {
    lm.out <- lm_robust(formula,
                          clusters = CountryName, data = FGMImpute$imputations[[i]])
    b.out <- rbind(b.out, summary(lm.out)$coefficients[,1])
    se.out <- rbind(se.out, summary(lm.out)$coefficients[,2])
    }
  combined.results <- mi.meld(q = b.out, se = se.out)
  mimp4.out <- combined.results

  mimp4.out <- data.frame(t(mimp4.out$q.mi), t(mimp4.out$se.mi))
  colnames(mimp4.out) <- c("coefs", "se")
  
  #Calculate test statistics
  mimp4.out$zscores <- mimp4.out$coefs/mimp4.out$se
  mimp4.out$zscores <- abs(mimp4.out$zscores)
  mimp4.out$pvalue <- 2*(1-pnorm(mimp4.out$zscores))
  return(mimp4.out)
  }

mi1 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + GDPpc + as.factor(CountryName))
mi2 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + popdensity + as.factor(CountryName))
mi3 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + netodagdp + as.factor(CountryName))
mi4 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + Trade + as.factor(CountryName))
mi5 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + enrollmentgenderparityindex + as.factor(CountryName))
mi6 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + autocracy + as.factor(CountryName))
mi7 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + durable + as.factor(CountryName))
mi8 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + instab + as.factor(CountryName))
mi9 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + war + as.factor(CountryName))
mi10 <- miresults(BirthYearAgeFGMCwtd ~ fgmclag + outlaw + as.factor(CountryName))

#Summary Statistics

library(fBasics)

stargazer(t(basicStats(fgmc)), out = "summary.html")